//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : subs_kernels.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef SUBS_KERNELS_H
#define	SUBS_KERNELS_H

namespace rrbsolver
{

template <class T>
__global__ void kernel_subs1(T* x, const Grid g)
{
    int ld = g.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;
    int u_r2 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_r2 = BORDER_WIDTH2 + g.cy + by * By + ty;
    int u_b1 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_b1 = BORDER_WIDTH + by * By + ty;
    int u_b2 = BORDER_WIDTH + bx * Bx + tx;
    int v_b2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_r1 = ld * v_r1 + u_r1;
    int loc_r2 = ld * v_r2 + u_r2;

    const T x_b1    = __FETCH_XX(u_b1,   v_b1,   ld);
    const T x_b2    = __FETCH_XX(u_b2,   v_b2,   ld);
    const T x_b1um1 = __FETCH_XX(u_b1-1, v_b1,   ld);
    const T x_b2up1 = __FETCH_XX(u_b2+1, v_b2,   ld);
    const T x_b1vp1 = __FETCH_XX(u_b1,   v_b1+1, ld);
    const T x_b2vm1 = __FETCH_XX(u_b2,   v_b2-1, ld);

    const T c_b2    = __FETCH_CC(u_b2,   v_b2,   ld);
    const T c_b1    = __FETCH_CC(u_b1,   v_b1,   ld);
    const T c_b1um1 = __FETCH_CC(u_b1-1, v_b1,   ld);
    const T c_b2up1 = __FETCH_CC(u_b2+1, v_b2,   ld);
    const T c_b2vm1 = __FETCH_CC(u_b2,   v_b2-1, ld);;
    const T c_b1vp1 = __FETCH_CC(u_b1,   v_b1+1, ld);

    const T cx1 = c_b1 * x_b1;
    const T cx2 = c_b2 * x_b2;

    x[loc_r1] -= g.nn[loc_r1] * cx2 +                // north
                 g.ee[loc_r1] * cx1 +                // east
                 g.ss[loc_r1] * c_b2vm1 * x_b2vm1 +  // south
                 g.ww[loc_r1] * c_b1um1 * x_b1um1;   // west

    x[loc_r2] -= g.nn[loc_r2] * c_b1vp1 * x_b1vp1 +  // north
                 g.ee[loc_r2] * c_b2up1 * x_b2up1 +  // east
                 g.ss[loc_r2] * cx1 +                // south
                 g.ww[loc_r2] * cx2;                 // west
}

template <class T>
__global__ void kernel_subs2(T* x, const T* r, const Grid g)
{
    int ld = g.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;
    int u_r2 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_r2 = BORDER_WIDTH2 + g.cy + by * By + ty;
    int u_b1 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_b1 = BORDER_WIDTH + by * By + ty;
    int u_b2 = BORDER_WIDTH + bx * Bx + tx;
    int v_b2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_r1 = ld * v_r1 + u_r1;
    int loc_r2 = ld * v_r2 + u_r2;
    int loc_b1 = ld * v_b1 + u_b1;
    int loc_b2 = ld * v_b2 + u_b2;

    const T x_r2    = __FETCH_XX(u_r2,   v_r2,   ld);
    const T x_r1up1 = __FETCH_XX(u_r1+1, v_r1,   ld);
    const T x_r1    = __FETCH_XX(u_r1,   v_r1,   ld);
    const T x_r2um1 = __FETCH_XX(u_r2-1, v_r2,   ld);
    const T x_r2vm1 = __FETCH_XX(u_r2,   v_r2-1, ld);
    const T x_r1vp1 = __FETCH_XX(u_r1,   v_r1+1, ld);

    T sum;

    sum = r[loc_b1] - g.ss[loc_r2] * x_r2
                    - g.ee[loc_b1] * x_r1up1
                    - g.ss[loc_b1] * x_r2vm1
                    - g.ee[loc_r1] * x_r1;
    x[loc_b1] = sum * g.cc[loc_b1];

    sum = r[loc_b2] - g.nn[loc_b2] * x_r1vp1
                    - g.ww[loc_r2] * x_r2
                    - g.nn[loc_r1] * x_r1
                    - g.ww[loc_b2] * x_r2um1;
    x[loc_b2] = sum *g.cc[loc_b2];
}

} // namespace rrbsolver

#endif	/* SUBS_KERNELS_H */

